clc,clear;
close all;
%导入数据
data =xlsread("random_cities_coordinates.xlsx");
% 城市坐标
cities=data;
% 城市数量
nCities = length(cities);
[bP,bl,l]=LFACOs(nCities,cities);
disp('Best path:');
disp(bP);
disp('Best length:');
disp(bl);
maxIter = 100;          % 最大迭代次数
tusi(l,nCities,cities,bP,bl,maxIter);
%% 作图

function ktu=tusi(l,nCities,cities,bP,bl,maxIter)
l=[l(2:end)];
n=length(l);
count=1:n;
x=zeros(nCities+1,1);
y=zeros(nCities+1,1);
for i=1:nCities
    x(i)=[cities(bP(i),1)];
    y(i)=[cities(bP(i),2)];
end
x(i+1)=cities(bP(1),1);
y(i+1)=cities(bP(1),2);

figure
plot(count,l);
xlabel('迭代次数');
ylabel('最优值');
title('模拟迭代图');
hold on
figure;
plot(x,y,'-ro', 'LineWidth', 2, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'g', 'MarkerSize', 10);
xlabel('X 坐标')
ylabel('Y 坐标')
title(['Iteration', num2str(maxIter), ' - Best Length: ', num2str(bl)])
grid on; % 打开网格
end

%% ACO(蚁群算法)
function [bestPath,bestLength,L]= LFACOs(nCities,cities)

% 距离矩阵
dist = zeros(nCities, nCities);
for i = 1:nCities
    for j = 1:nCities
        if i ~= j
            dist(i, j) = norm(cities(i, :) - cities(j, :));
        else
            dist(i, j) = inf; % 自己到自己的距离设为无穷大
        end
    end
end

% 参数设置
nAnts = 50;             % 蚂蚁数量
maxIter = 100;          % 最大迭代次数
alpha = 1;              % 信息素重要性
beta = 5;               % 距离重要性
rho = 0.9999;           % 信息素挥发系数
Q = 10;                 % 信息素增加常数
alapa=1.8;              %莱维飞行幂律分布的指数(0~2)
albt=0.1;               %初始步长
% 初始化信息素矩阵
pheromone = ones(nCities, nCities);
% 记录最佳路径及其长度
bestPath = [];
bestLength = inf;
L=[];
% figure;
% 迭代优化
for iter = 1:maxIter
    % 每只蚂蚁选择路径
    paths = zeros(nAnts, nCities);
    lengths = zeros(nAnts, 1);
    for k = 1:nAnts
        paths(k, :) = conS(pheromone, dist, alpha, beta, nCities);
        lengths(k) = calul(paths(k, :), dist);
    end

    % 更新最佳路径
    [minLength, idx] = min(lengths);
    L=[L,bestLength];
    
    if minLength < bestLength
        bestLength = minLength;
        bestPath = paths(idx, :);
    end
    if rho<=0||rho>=1
        rho=0.6;
    end
    % 更新信息素
    pheromone = (1 - rho) * pheromone;
    for k = 1:nAnts
        for i = 1:nCities-1
            pheromone(paths(k, i), paths(k, i+1)) = pheromone(paths(k, i), paths(k, i+1)) + Q / lengths(k);
        end
        pheromone(paths(k, nCities), paths(k, 1)) = pheromone(paths(k, nCities), paths(k, 1)) + Q / lengths(k);
    end
    % 绘制当前最优路径
%     plot(cities(best_path, 1), cities(best_path, 2), 'o-');
%     title(['Iteration ', num2str(iter), ' - Best Length: ', num2str(bestLength)]);
%     pause(0.01); % 暂停一段时间以便观察图形变化
    %引用莱维飞行改变rho从而提高搜索能力
    Lr = levyflight(alapa);%分布步长
    rho=rho-Lr*albt;
%     albt=albt*(1/iter)^0.1;
end
    bestPath=[bestPath,bestPath(1)];%回到起点
end
%% 选择路径
function path = conS(pheromone, dist, alpha, beta, nCities)
    path = zeros(1, nCities);
    visited = false(1, nCities);
    path(1) = 5; % 选择起始城市
    visited(path(1)) = true;
    for i = 2:nCities
        currentCity = path(i-1);
        probabilities = calP(currentCity, visited, pheromone, dist, alpha, beta, nCities);
        path(i) = rouS(probabilities);
        visited(path(i)) = true;
    end
end
%%   信息素计算
function probabilities = calP(currentCity, visited, pheromone, dist, alpha, beta, nCities)
    probabilities = zeros(1, nCities);
    for j = 1:nCities
        if ~visited(j)
            probabilities(j) = (pheromone(currentCity, j)^alpha) * ((1 / dist(currentCity, j))^beta);
        end
    end
    probabilities = probabilities / sum(probabilities);
end
%%  选择
function selected = rouS(probabilities)
    cumulativeSum = cumsum(probabilities);
    r = rand();
    selected = find(cumulativeSum >= r, 1);
end
%%  目标值计算
function length0 = calul(path, dist)
    length0 = 0;
    nCities = length(path);
    for i = 1:nCities-1
        length0 = length0 + dist(path(i), path(i+1));
    end
    length0 = length0 + dist(path(nCities), path(1)); % 回到起点
end

% 莱维飞行生成函数
function L = levyflight(beta)
    % 参数
    sigma_u = abs((gamma(1 + beta) * sin(pi * beta / 2) / (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2))))^(1 / beta);
    sigma_v = 0.5;

    % 生成莱维分布步长
    u = randn * sigma_u;
    v = randn * sigma_v;
    L = u / abs(v)^(1 / beta);
end
